Tagged-particle dynamics in a hard-sphere system: mode-coupHng theory analysis 
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The predictions of the mode-coupling theory of the glass transition (MCT) for the tagged-particle 
density-correlation functions and the mean-squared displacement curves are compared quantitatively 
and in detail to results from Newtonian- and Brownian-dynamics simulations of a polydisperse quasi- 
hard-sphere system close to the glass transition. After correcting for a 17% error in the dynamical 
length scale and for a smaller error in the transition density, good agreement is found over a wide 
range of wave numbers and up to five orders of magnitude in time. Deviations are found at the 
highest densities studied, and for small wave vectors and the mean-squared displacement. Possible 
error sources not related to MCT are discussed in detail, thereby identifying more clearly the issues 
arising from the MCT approximation itself. The range of applicability of MCT for the different types 
of short-time dynamics is established through asymptotic analyses of the relaxation curves, exam- 
ining the wave-number and density-dependent characteristic parameters. Approximations made in 
the description of the equilibrium static structure are shown to have a remarkable effect on the 
predicted numerical value for the glass-transition density. Effects of small polydispersity are also 
investigated, and shown to be negligible. 



PACS numbers: 64.70.Pf, 82.70.Dd 



I. INTRODUCTION 



Understanding the slow dynamical processes that oc- 
cur when one cools or compresses a liquid is a great 
challenge of condensed matter physics. In particular in 
the time window accessible to scattering experiments or 
molecular-dynamics (MD) computer simulations, one ob- 
serves in equilibrium a precursor of the liquid-glass tran- 
sition that is commonly termed structural relaxation. 
From these experiments, a large amount of detailed in- 
formation about the equilibrium fluctuations in such sys- 
tems is available [ij. 

Many of the recent experiments on structural relax- 
ation were stimulated by the mode-coupling theory of the 
glass transition (MCT). This theory attempts to provide 
a first-principles description of the slow structural relax- 
ation processes, requiring as input the (averaged) equi- 
librium static structure of the system under study. Un- 
fortunately, for many commonly studied glass formers, 
the latter is not available to the extent required. Thus 
comparisons of MCT with experiment usually proceed 
by referring to asymptotic predictions or schematic sim- 
plifications of the theory that can be evaluated without 
restriction to a specific system, and by fitting the re- 
maining parameters of the theory. One has to be careful 
when interpreting these results, since it is known that 
most experimental data is hardly inside the regime of 
applicability of the asymptotic formulas [3,l3l- Still, this 
way, many studies of the predicted MCT scenario have 



been performed (see Ref. [l| for a review). 

Having established the general scenario, important 
questions arising are what are its ranges of validity, and 
what is the effect of the approximations made in the 
course of deriving the theory. These questions can be 
addressed by comparing the 'full' solutions of the theory 
to experimental results for one and the same system for 
which the static structure is known in detail. While work 
has been done along this direction to various degrees of 
detail recently, a coherent picture for a single prototyp- 
ical system has not yet emerged. This paper aims to- 
wards filling this gap by providing a detailed comparison 
of computer-simulation results for a system of quasi-hard 
spheres with the corresponding 'full' MCT solutions, to 
establish the performance of MCT in describing the dy- 
namics of a prototypical glass-forming system as a fully 
microscopic theory. 

Such first-principles comparisons have become possible 
with the appearance of powerful MD simulations for sim- 
ple model systems. Simulation data has been used to suc- 
cessfully test the MCT predictions for the frozen glassy 
structure (the long-time limit of the dynamical correla- 
tion functions) for a mixture of Lennard-Jones particles 
U, a model silica melt ^, and a computer model of 
the molecular glass former Ortho-terphenyl Jjj. In these 
cases, the equilibrium-structure input to MCT was de- 
termined from the simulations themselves. The dynami- 
cal information has not been compared to MCT in these 
cases. This comparison has been tackled for the Lennard- 
Jones mixture [3 and for two binary hard-sphere mix- 



tures J3|, but there the discussion had to be restricted 
to the slowest decay process, while qualitative deviations 
from MCT at intermediate and short times could not 
be resolved. This is in contrast to a fuU-MCT analysis 
of experimental light-scattering data from a quasi-binary 
hard-sphere like colloidal mixture |2|, where agreement 
over the full accessible time window was found as far as 
MCT was concerned, including short and intermediate 
times. It is unclear to what extent the different system 
types and the different forms of short-time dynamics be- 
tween the MD simulations and the colloidal system give 
rise to the differing results. Thus it seems appropriate to 
perform this comparison for an even more fundamental, 
paradigmatic glass former, viz.: the hard-sphere system 
(HSS). 

Simulations for this system close to the glass transi- 
tion have been performed by Doliwa and Heuer 10, 11] 
using a Monte Carlo procedure and a slight polydis- 
persity There, an emphasis was put on the analysis 
of cooperative motion on the single-particle level, and 
no quantitative connection to MCT was reported. In- 
stead we focus on the analysis of the self-intermediate 
scattering functions, which can be directly compared to 
theory and experiment. We chose to perform molecu- 
lar dynamics (MD) simulations instead of MC, in order 
to be able to also study the influence of different realis- 
tic types of short-time dynamics, i.e. 'atomistic' Newto- 
nian dynamics (ND) and 'colloidal' Brownian dynamics 
(BD). Such a study has been performed earlier for the 
Lennard- Jones mixture mentioned above |12| . however 
no fuU-MCT analysis was included there. 

For an observation of the equilibrium glassy dynam- 
ics, it is, in general, necessary to avoid crystallization by 
some means. For the HSS, this can be accomplished by 
introducing a small amount of polydispersity that drasti- 
cally reduces crystallization rates 13] . This is inherently 
the case in studies of colloidal suspensions. In the MD 
simulation, we are able to fully control the distribution 
of particle radii in the system. In colloidal suspensions, 
solvent-mediated hydrodynamic interactions (HI) are in- 
evitable. It is an as yet not settled question to what 
extent HI influence the dynamics at high densities. In 
the present simulations, HI are not present. Thus our 
study also serves to complement previous analyses of col- 
loidal hard-sphere suspensions with asymptotic formulas 
|14L Il5| . demonstrating that HI are not an important 
ingredient for a quantitative description of structural re- 
laxation. 

The paper is organized as follows. First, we introduce 
in Sec. ^ the relevant quantities for the discussion. An 
investigation of some asymptotic properties of the sim- 
ulation data is performed in Sec. IIIII whereas Sec. IIVI 
is devoted to a comparison of the time-dependent data 
with MCT results for the one-component HSS. In Sec.lVl 
the effects of polydispersity will be discussed within the 
framework of this MCT analysis. We summarize our find- 
ings in the conclusions. Sec. IVII 



II. SIMULATION AND THEORY DETAILS 

A. Molecular Dynamics Simulations 

We perform standard molecular-dynamics simulations 
of A^ = 1000 particles in the canonical ensemble in a 
polydisperse system of quasi-hard spheres. The core-core 
repulsion between particles at a distance r is given by 



14 (r) = /cbT — 



-36 



(1) 



where di2 is the center-to-center distance, di2 — {di + 
^2)72, with di and ^2 the diameters of the particles. This 
potential is tailored to be a continuous approximation 
to the hard-sphere potential considered in the theoret- 
ical part of the work, as this facilitates the simulation 
of Brownian dynamics. The control parameters of this 
soft-sphere system are the number density g and temper- 
ature T; they appear however only as a single effective 
coupling parameter, F = gT~^^ 16]. In the simulation, 
F is varied by changing the density and keeping the tem- 
perature fixed. In the following, we denote the number 
density in terms of a packing (or volume) fraction, which 
for a monodisperse system reads (p = {TT/6)gd'^. To avoid 
crystallization, the diameters of the particles in the sim- 
ulation are distributed according to a flat distribution 
centered around d and a half- width of 6/2 — O.ld. Thus 
the volume fraction reads p — (7r/6)(i^ [l + S^] g. 

Note that, due to polydispersity and finite-size effects, 
it is not trivial to ensure that the volume fraction remains 
constant among different runs, i.e. different realizations 
of the polydispersity distribution. If one randomly draws 
N particles with radii according to the polydispersity dis- 
tribution at a fixed number density, the resulting packing 
fraction will vary from run to run, by up to about 1% in 
the cases we have investigated. This is not acceptable, 
since the slow dynamics to be discussed depends sensi- 
tively on the packing fraction. In order to eliminate such 
fluctuations of ip, we instead choose a fixed realization of 
the radius distribution (1000 equally spaced radii from 
0.9 to 1.1), and randomly assign each radius to one of 
the particles in the initial configuration. 

Both Newtonian and Brownian dynamics simulations 
were performed, to analyze the effect of the microscopic 
dynamics on the structural relaxation. Newtonian dy- 
namics (ND) was simulated by integrating the Newton's 
equations of motion in the canonical ensemble at constant 
volume. In Brownian dynamics (BD), or more precisely, 
strongly damped Newtonian dynamics, each particle ex- 
periences a Gaussian distributed white noise force with 
zero mean, ff(t) , and a damping force proportional to the 
velocity, "ff, apart from the deterministic forces from the 
interactions. Hence the equation of motion for particle j 



-E^. = - 



jrj+rfj{t), 



(2) 



where 7 is a damping constant. The stochastic and 
friction forces fulfill the fluctuation-dissipation theorem, 
{f]iit)f]jit')) = 6kBT"fS(t - t')S^j. The value of 7 was 
set to {30/^/3)kT/{dvth)', this 'overdamped limit' en- 
sures that the results presented here no longer show a 
dependence on the value 7. Such a form of the dynam- 
ics has been introduced in the study of glassy relaxation 
by Gleim et al. ^3|■ Let us note that the short-time 
dynamics visible in the correlators and in the mean- 
squared displacement still is not strictly diffusive, but 
rather strongly damped ballistic. Since it is not our aim 
to investigate the very short-time dynamical features of 
the simulations, this will not be discussed in the follow- 
ing. 

Equilibration runs were done with ND in all cases, 
since the damping not only introduces a change in the 
overall time scale but also slows down the equilibration 
process. Lengths are measured in units of the diameter, 
the unit of time is fixed setting the thermal velocity to 
I'th = \/k-BT/m = l/v3, and the temperature is fixed 
to A;bT = 1/3. In ND, the equations of motion were in- 
tegrated using the velocity- Verlet algorithm 17], with a 
time step St — 0.0025. The thermostat was applied by 
rescaling the particle velocity to ensure constant temper- 
ature every nt time steps. For well equilibrated samples, 
no effect of nt was detected. The equations for Brownian 
motion were integrated following a Heun algorithm [13 
with a time step of St = 0.0005. In this case, no exter- 
nal thermostat was used, since the samples were already 
equilibrated when running BD simulations. 



The orientational order parameter Qq |19ll2n| was used 
to check that the system was not crystalline. For amor- 
phous liquid-like structures, Qq is close to zero, whereas 
it takes a finite value for an ordered phase. Even though 
the 10% polydispersity is high enough to prevent crystal- 
lization in most cases, some samples at the highest vol- 
ume fractions studied did show a tendency to crystallize. 
Those have been excluded from the analysis. 

The volume fractions investigated in this work are 
ip = 0.50, 0.53, 0.55, 0.57, 0.58, 0.585, and 0.59. At 
each volume fraction, we extracted as statistical in- 
formation on the slow dynamics the self part of the 
intermediate scattering function for several wave vec- 
tors (f, (/)'*(g,i) — (exp[— i(f(rs(i) — rs(0))]), formed with 
the Fourier-transformed fluctuating density of a single 
'tagged' particle at position f\j(i). Here, angular brack- 
ets (•) denote canonical averaging. A related quan- 
tity which we also extracted from the simulations is the 
mean-squared displacement (MSD) of a tagged particle, 
Sr'^it) = {{rs{t) - fs(0))^). The correlators and the MSD 
were averaged over typically 50 runs, except for the BD 
simulations at 1^ = 0.585 and 0.59, where 20 runs have 
been performed originally. For ip = 0.59 we have also 
performed additional runs for both ND and BD in order 
to investigate some phenomena found there, see Sec. lIIIBl 
below. 



B. Mode-Coupling Theory 

In a system of N structureless classical particles, i.e. 
without any internal degrees of freedom, the statistical 
information on the structural dynamics is encoded in 
the density correlation function, (/)(g,t) — {Q{q,t)*g{q)), 
formed with the fluctuating number densities g{q, t) = 
^- -^exp(ig • fj{t))/^/N for wave vector q. 4>{q,t) is a 
real function that depends on g only through q — \q\, 
since it is the Fourier transform of a real, translational- 
invariant and isotropic function. The dynamics given 
by (j){q,t) is probed by the mean-squared displacement, 
Sr'^{t), and the self-part of the intermediate scatter- 
ing function (also called tagged-particle correlation func- 
tion), (/'"(g, i), extracted from our simulations. Note that 
the latter is linked to the MSD in the limit (7 — > 0, via 
<j>^iq,t) = l-{l/6)q^Sr^t) + 0{q^). 

The mode-coupling theory of the glass transition 
(MCT) 1^ builds upon an exact equation of motion for 
the density autocorrelation function (l>[q,t), 



-J-^d^^iq,t)+<j)iq,t) 

+ m{q,t~t')dt'(t){q,t')dt' = 0. (3a) 
Jo 

Here, ^{q)^ — q^v1^/S{q) is a characteristic squared fre- 
quency of the short-time motion. The equation of motion 
is supplemented by the initial conditions (f){q, i = 0) = 1 
and dt4>{q, i = 0) = 0. All many-body interaction effects 
are contained in the memory kernel m{q,t), the descrip- 
tion of which is the aim of the MCT approximations. One 
splits off from this kernel a mode-coupling contribution 
rn?^'~^'^{q,t), while the remainder is assumed to describe 
regular liquid-state dynamics. Let us approximate this 
latter part by an instantaneous contribution. 



m{q,t)^ii,{q)/niqy)5{t)+m^'^'{q,t) 



(3b) 



The damping term ^{q) is chosen as i^ = (30/\/3) vth/d 
independent of g; a choice that ensures the short-time 
expansion of (f){q, t) in the overdamped limit to match 
that one of the simulation, cf. Eq. Q: one gets (j){q, t) = 
1 - {qyS{q)){kBT/j)t + 0{t^) = 1 - {n{qf/v)t + 0{t^). 
Note that the g-independent choice of v destroys mo- 
mentum conservation for the hard-sphere particles; this 
is appropriate for a model of a colloidal system. 

The MCT contribution to the memory kernel is given 
by m^^'^{q,t) = Tq[(j){t)], where 



m 



Q 
2g4 



—-S{q)S{k)S{p)V{q,k,p)f{k)f{p) 



(3c) 

and the abbreviation p — q — k has been used. The ver- 
tices V(q, k,p) are the coupling constants of the theory, 
through which all crucial control-parameter dependence 
enters. They are given entirely in terms of static two- and 



three-point correlation functions describing the equihb- 
rium structure of the system's hquid state. The latter 
are approximated using a convolution approximation, so 
that the vertex reads 



V{q, k,p) = {qk)c{k) + {qp)c{p) 



(3d) 



Here, c{q) is the direct correlation function (DCF) con- 
nected to the static structure factor by S{q) = 1/(1 — 
gc{q)). 

The long-time limit of the correlation functions, f{q) = 
limt_>oo (t>iq^ t), is used to discriminate between liquid and 
glassy states. In the liquid, f{q) = 0, while the glass is 
characterized by some f{q) ^ 0. From Eqs. lO, one 
finds f{q) as the largest real and positive solution of the 
implicit equations i22j 



/(g) 



•^J/] 



(4) 



In particular, there exist critical points in the control- 
parameter space, identified as ideal glass transition 
points, where a new permissible solution of Eq. Q ap- 
pears. Typically, f{q) jumps discontinuously from zero 
to nonzero values there. Close to such a critical point on 
the liquid side, the correlation functions exhibit a two- 
step relaxation scenario, composed by a relaxation to- 
wards a plateau value, and by a later relaxation from this 
plateau value to zero termed a relaxation. On approach- 
ing the transition, the characteristic time scale for the 
a relaxation diverges, and an increasingly large window 
opens where the correlation functions stay close to their 
plateau. The plateau values on the liquid side are in lead- 
ing order given by the critical solutions of Eq. ^ , f^ {q) , 
i.e. by the maximal solutions of Eq. Q) evaluated at a 
critical point. The time window for which 4>{q,t) is close 
to f^iq) is called the /3-relaxation regime, and is the ob- 
ject of asymptotic predictions of MCT HHlllil. These 
include scaling laws for the correlators, whose power-law 
exponents a and 6, called the critical and the von Schwei- 
dler exponent, are given by an exponent parameter A. 
The latter is calculated within MCT and depends on the 
static equilibrium structure of the system. We will test 
some of the predictions connected with (3 relaxation in 
Sec lilTDl 

Let us also recollect the MCT equations of motion 
for the tagged-particle correlation function (j)''{q,t) of a 
tagged particle that is of the same species as the host 
fluid, since this will be the quantity we shall analyze be- 
low. For it, an expression similar to that of Eqs. © 
holds, 

:dfct>'{q,t) + (t>'{q,t) 



+ f m''{q,t-t')df<j)'{q,t')dt' ^Q, (5a) 
Jo 

where we have ri*(q)^ = q^v^^. The tagged-particle mem- 
ory kernel is given in MCT approximation by m''{q, t) « 



iUsiq)/Omega-'{qf)5it) + T'[(j)'{t),cj){t)], with 
1 f fPk 

•^1/^/] = 71 1 77^v%q,k)f{k)r{p) , (5b) 



q'^ J (27r)3 



and with vertices 



V'{q,k)^iqkfc{ky^ 



(5c) 



where we set Vsiq) = J^ in the following. The quali- 
tative features of (j)''{q,t) close to an ideal glass transi- 
tion are the same as those of (j){q,t), as long as it cou- 
ples strongly enough to via Eq. I)5b|l. In this generic 
case, also 4>^{q, t) develops a two-step relaxation pattern, 
with plateaus given by the critical solution f^''^{q) of the 
tagged-particle analog of Eq. Q , 



r(g) 

i-/-^(g) 



= K[.f,h^ 



(6) 



The mean-squared displacement (MSD) 6r'^{t) can be 
calculated from the g — > limit of the tagged-particle 
correlation function. One gets 

dt5r\t) + vl^ [ mg(i - t')Sr^t') dt' = &vl^t , (7) 
Jo 

where we have set TOg(i) = limg^o q'^rn^{q,t). 

Eqs. (Q can be solved numerically for the functions 
(f){q,t), once the vertices V{q,k,p) have been calculated 
from liquid-state theory. To this end, the wave vectors 
are discretized on a regular grid of M wave numbers with 
spacing A^: qd — lAq + qo- We have used M = 300, 
Aq = 0.4/3, and go — 0.2/3, implying a cutoff wave vec- 
tor q*d = 39.93. This discretization is enough to ensure 
that the long-time part of the dynamics does not show 
significant numerical artifacts '^ and has been used be- 
fore in the discussion of MCT results for the HSS [2J|. 
Once the (j){q, t) have been determined, a similar numer- 
ical scheme allows to evaluate Eqs. for the 4)^{q,t), 
and from this, one gets 5r^{t) from Eq. {Tj). 

For the solution of Eqs. Q, a straightforward itera- 
tion scheme guarantees a numerically stable determina- 
tion of the correct solutions /(g) = (/>(g, t -^ oo) 22] and, 
once the f{q) are calculated, of f'iq) = (fi^iq^t — > oo). 
From the distinction between states with f{q) ^ and 
/(g) = 0, the critical point •-p'^ can be found by iteration 
in iy9. For the solution of these equations, we have used 
a discretization with M = 100, Ag = 0.4, and go = 0.2. 
This is sufficient to ensure that errors in the /(g), /"(g), 
and Lp'^ resulting from the different discretizations used 
are small. 

A few results shall also be discussed concerning the 
polydispersity- induced effects. MCT for continuous poly- 
dispersity distributions is not available, but we try to es- 
timate the influence of the polydispersity by calculating 
MCT results for S'-component mixtures with the species' 
diameters chosen to mimic the simulated polydisperse 
distribution. We have used an 5* = 3 model with diame- 
ters da e {1 — w, 1, 1 -I- ui}, and Qa — £'/3, where a labels 



the species of the mixture, and Qa is the partial number 
density of each species. Here, we set w — l/\/200 in 
order to match the second moment of the discrete dis- 
tribution to that of the one used in the simulation. We 
have also calculated results for an S" = 5 model, with 
da e {0.9,0.95,1.0,1.05,1.1}, and g^ = g/5, chosen to 
contain particles within the same size range as in the 
simulation. The MCT equations, Eqs. Q, generalize to 
mixtures in an obvious way, leading to equations of mo- 
tion for the matrix of partial density correlators, ^affiq, t) 
[23,123. Similar to Eqs. Q, the correlators for a tagged 
particle of either one of the species, (/)^(q,i), are calcu- 
lated, together with their long-time hmits, /„('?)■ We can 
now define 'averaged' tagged-particle quantities as 



/pd(9) = ^E/'^(9)' 



(8) 



a=l 



and similarly for (f)t^^{q^t). These quantities are analo- 
gous to the quantities extracted from the polydisperse 
MD simulations. 

To calculate results from the MCT equations, we re- 
quire as the only input expressions for the direct corre- 
lation function c{q) entering the vertices, Eqs. Ij3d|l and 
(|5c|. For the multi-component analog of these expres- 
sions, one requires knowledge of the full matrix of direct 
correlation functions, Cai3{q). The DCF could be either 
determined from simulations, or taken from well-known 
results of liquid structure theory. For hard spheres, 
the Percus-Yevick (PY) closure to the Ornstein-Zernike 
equation pro vides a fairly accurate parameter-free de- 
scription ^^. Using the PY-DCF as input to MCT, we 
thus obtain results for the dynamics of the HSS that are 
independent on any empirical parameters or any (usually 
not readily available) simulation input. These predic- 
tions are the best we can currently achieve from within 
MCT as a completely parameter-free theory. Further- 
more, the wealth of asymptotic predictions of MCT has 
been worked out in great detail for this model [j, [23 ■ 
Note that the PY approximation to the DCF itself intro- 
duces errors that are independent from those introduced 
by the MCT approximation. It has been pointed out 
recently that these PY-induced errors can be quite pro- 
nounced in the MCT-calculated quantities, even if they 
appear small at the S{q) level Q- To disentangle these 
two error sources, we have also performed some calcula- 
tions within MCT with S{q) obtained from our simula- 
tion, as will be discussed below. 



III. DATA ANALYSIS 

Let us start the discussion of the data by a comparison 
of the structure factor S{q) obtained from the simulation 
with the PY approximation, since this is the crucial input 
to all MCT calculations below. Fig.Qshows this quantity 
for (fi = 0.58. While PY reproduces the oscillation period 
in S{q), i.e. the typical length scale, correctly, it overes- 




FIG. 1: Comparison of the static structure factor S{q) for 
the simulated polydisperse soft-sphere system at <^ = 0.58 
(symbols) with the Percus-Yevick approximation to the hard- 
sphere S{q) at the same density (solid line). The dashed line 
shows the Percus-Yevick result for ip = 0.505. The region 
around the first diffraction peak is enlarged in the inset. 
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FIG. 2: MCT results for the critical nonergodicity parame- 
ters f^iq), using as input the static structure factor S{q) from 
Percus-Yevick theory for hard spheres (dashed line). The 
crosses connected by the solid line show the results using S{q) 
as obtained from the polydisperse nearly-hard-sphere simula- 
tion. 



timates the peak heights, i.e. the strength of ordering in 
the system llq| . Since the strength of the coupling con- 
stants in MCT is directly connected to the peak heights 
in S{q), the MCT calculation based on the PY S{q) will 
overestimate the tendency to glass formation. One can 
try to adjust the peak heights by setting a lower pack- 
ing fraction in the PY calculation. This is demonstrated 
by the dashed line in Fig. ^ where (p = 0.505 has been 
taken. This value has in fact been determined by the 
MCT fits presented in Sec. IIVI and is chosen such that 
the final relaxation time in the MCT calculations at that 
density matches the one of the simulations at (/? = 0.58. 
As Fig. n demonstrates, this introduces a small error in 
the oscillation period in S{q). 

It is well known that MCT, based on the PY struc- 
ture factor for hard spheres, underestimates the glass- 
transition packing fraction of that system. One gets 
V'mct ~ 0.516 13, instead of the value reported from 
experiments on colloidal hard-sphere-like suspensions. 



ip'^ K, 0.58 [lj|. In order to determine to which extent 
such an underestimation can stem from deviations of PY 
from the simulated S'(q), which are visible in Fig. ^ we 
have calculated MCT results for lySMCT ^'^^ the critical 
plateau values /"^(g) both using the PY approximation 
and using our simulation results for S'(g) as input to the 
theory. We have evaluated the structure factor from the 
simulation at (/? = 0.50 and ^p = 0.58, where we could 
get reasonable statistics for this quantity. The MCT cal- 
culations then proceed by a linear interpolation between 
these two cases to approximate S{q) at nearby values of 
ip. The critical nonergodicity parameters f'^{q) thus ob- 
tained are shown in Fig. [5] They agree well for qd > 6, 
lending confidence to the PY-based discussion of the cor- 
relation functions. Smaller q have been omitted from the 
figure since there, differences can be seen that are related 
to insufficient statistics for the simulated S{q) at small 
wave numbers. The results for the exponent parameter 
A also do not differ significantly between the two cal- 
culations. We get A ~ 0.735 in the PY-based case 0, 
and 0.727 < A < 0.773 based on the simulated S{q), 
the latter value depending somewhat on the discretiza- 
tion used. The values for the critical packing fraction, 
however, differ between the two calculations: instead of 
'/'mct ~ 0.516, we get </?mct ~ 0.585 when using the 
simulated data to obtain S{q). Note that this makes 
this MCT result almost coincide with what has been re- 
ported for colloidal realizations of a hard-sphere system 
[Tjl- Such agreement is accidental, particularly because 
the value 'p'^im. extracted from our simulations is even 
higher, but it demonstrates that the approximations used 
for S{q) need critical assessment. Let us also note that 
the findings described here do not completely agree with 
similar results reported in Ref. j^]- There, the same qual- 
itative trend for (p^cx ^^^ found for a hard-sphere sys- 
tem, and as well for two binary hard-sphere mixtures. 
But in this case, the values for /'^(q) based on the sim- 
ulated structure-factor input differed notably from those 
calculated within the PY approximation, while we find 
no significant difference in this quantity. In principle, 
our simulation-based results for f'^{q) have no reason to 
be closer to the PY results than the simulation-based re- 
sults from Ref. [g, since we use a slightly polydisperse 
soft-sphere system, while in Ref. Q, strictly monodis- 
perse hard spheres have been simulated, at the cost of 
having to extrapolate to the desired high densities. 

The results shown in Fig. |21 suggest that we may pro- 
ceed in the following discussion by basing all MCT results 
on the Percus-Yevick approximation for S{q). While this 
will make an adjustment of packing fractions (^mct nec- 
essary, it has the advantage of giving a first-principles 
theory to compare the simulation data to. In particular, 
the results presented above point out that neither the 
shape and strength of the a relaxation, nor the asymp- 
totic shape of the correlators in the /3-scaling regime will 
change much between the PY-based results and those 
based on the simulated structure factor. 

Before we embark on the comparison of the interme- 




FIG. 3: Simulation results for tp = 0.50, 0.53, 0.55, 0.57, 0.58, 
0.585, 0.59 (from left to right) at wave vector qd = 7.8. Heavy 
solid lines are the results using Brownian dynamics, thin lines 
the results for Newtonian dynamics. For the latter curves, 
times t have been multiplied by factors i* given in the inset. 
The dotted line is the BD result for tp = 0.01, indicating the 
dilute limit of the correlation function. The solid diamonds 
indicate the points where Newtonian and Brownian dynamics 
results start to agree at a 2% level. 



diate scattering functions with the 'full' MCT solutions, 
let us first analyze the simulation data according to the 
asymptotic predictions of MCT, in oder to identify the 
time window where MCT should certainly be applicable. 



A. Identification of Structural Dynamics 

In Fig. 131 results of the simulations are shown for dif- 
ferent packing fractions p). A wave vector qd = 7.8 close 
to the first peak in the static structure factor has been 
chosen. Different values of q show qualitatively similar 
scenarios. The thick solid lines in the figure are the simu- 
lation results for 'Brownian' dynamics simulations. Upon 
increasing (/3, one observes the emergence of a two-step 
relaxation process at times long compared with typical 
liquid time scales. A typical relaxation curve for the di- 
lute case, is exemplified by the dotted line in the figure, 
showing the BD simulation result for p) = 0.01. From 
this, we read off a 'microscopic' relaxation time for the 
short-time relaxation of t « 1. The slow two-step re- 
laxation pattern is usually referred to as structural re- 
laxation and is a precursor of the approach to a glass 
transition at some p'^. The scenario has been observed 
repeatedly in similar systems. In our simulations, we are 
able to follow the structural relaxation scenario for up to 
about five orders of increase in the relaxation time. 

MCT predicts that the structural relaxation becomes, 
up to a common time scale toj independent on the type of 
microscopic motion that governs the relaxation at short 
times. To demonstrate that this is the case. Fig. O 
shows as thin lines the simulation results using Newto- 
nian short-time dynamics. The data have been scaled in 
t in order to match the BD data at corresponding pack- 
ing fractions and at long times. Indeed, then the relax- 
ation curves match within our error bars at times t > 10, 





FIG. 4: Comparison of the BD data for (^ = 0.58 (plus 
symbols) at wave vectors qd = 4.0, 7.8, 13.8, and 19.8 (from 
top to bottom). The long-time part of the data for tp = 
0.55 (I symbols), 0.57 (crosses), 0.585 (circles), and 0.59 (star 
sybmols) is also shown, scaled in t to match the long-time 
part of the tp — 0.58 data. Solid lines are the MCT master 
curves for shifted qmct (see text for details). 



indicated by the diamond symbols in Fig. Only at 
shorter times, the regime of non-structural relaxation can 
be identified by the different shapes of the BD and ND 
curves. According to MCT, the scale factor t* = t^^ /t^^ 
used to match the BD and ND data at long times should 
be a smooth function of if, given by a constant in leading 
order close to f'^. For our simulation results, the values 
are as shown in the inset of Fig. O they are compatible 
with a constant shift t* ~ 4.25 within error bars. Only at 
ip = 0.58 and ip = 0.585 do we note a stronger deviation, 
the reason of which is unclear. The overall variance in 
t,t(</j) is comparable to the one found in a similar study 
of a binary Lennard- Jones mixture |l2l | . 

We conclude that the time window t > 10 deals 
with structural relaxation and thus comprises the regime 
where MCT predictions can be tested. At shorter times, 
deviations from those predictions must be expected. 
There are indications from theory [23, |23| and ND sim- 
ulations T| that these deviations are larger in ND. We 
thus primarily discuss the BD data, which prove to be 
simpler to understand within an MCT description. 



B. Qf-process analysis 

The second step of the two-step relaxation process, i.e. 
the decay of the correlators from their plateau value, is 
referred to as the a process. A prediction of MCT is 
that the shape of this a relaxation becomes independent 
on Lp in the limit Lp ^ p'^ — Q. Thus, scaling the corre- 
lation functions for given q and different pi to agree at 
long times should collapse the data onto a master curve. 
Figure 01 demonstrates the validity of a scaling for the 
BD data at several wave vectors between q — 4.0 and 
19.8. The scaling works as expected from the MCT dis- 
cussion of the HSS U for pi < 0.58. The closer a state 
is to (p'^, the larger is the window where the correlator 



FIG. 5; Demonstration of the variability between different 
simulation runs for the ND and BD simulations at p> = 0.59. 
Data is shown for qd = 7.8, with open (filled) symbols de- 
noting BD (ND) results. Triangles indicate averages over a 
small subset of the data (8 out of 30 sets for BD; 25 out of 
70 for ND) only; inverted triangles are the averages over the 
remaining data sets. The solid lines without symbols are the 
total averages. Dotted lines indicate the time-scaled MCT 
a-master functions. 



follows the a-master function. The increase of the (f>{q, t) 
above the master functions at shorter times is connected 
to the (3 process, discussed below. For (p = 0.59, a scal- 
ing breaks down at long times, t > 500. The reason for 
this is unclear, and cannot be understood within MCT. 
As observed by the orientational order parameter Qq, 
the system did not show appreciable trends to crystal- 
lization in any of the analyzed simulation runs. Also for 
ip — 0.585, some deviations from a scaling can be seen, 
particularly at g = 7.8 and at around t — 1000, which are 
not in agreement with the pre-asymptotic corrections to 
MCT a scaling. But in this case, the deviations are less 
pronounced than those at ip = 0.59. 

The behavior of the long-time dynamics at these two 
densities, p> = 0.585 and 0.59 is not fully understood. 
We have tried to improve the statistical averaging by in- 
creasing the number of simulation runs. However, there 
appear to be two subsets among the runs: one where the 
a-scaling violation is very pronounced, and one where 
the correlators instead follow the scaling behavior much 
closer. This happens in both the ND and the BD simula- 
tions, although the effect is more clearly seen in the ND 
case. Out of the 30 data sets we have averaged in the BD 
case for tp = 0.59, only 8 show the scaling violation; in the 
ND case we have averaged over 70 sets, with 25 of them 
deviating from scaling. While Fig. ^ shows the data av- 
eraged over all simulation runs. Fig. |S1 demonstrates the 
variation in a-time scale between the two types of data 
sets, obtained by restricting the averaging to the number 
of data sets specified above. While we have no a priori 
justification to modify the averaging procedure in any 
way, it allows us to point out, that possibly some 'rare' 
events take place in the system at these high densities, 
which we cannot classify as crystallization events on the 
basis of Qe, but which modify the dynamical long-time 
behavior in a distinct way. For the ND data, the a-time 



scale varies by a factor of 2.5 between the two cases. In 
the BD data, the effect is less pronounced, but still gives 
a factor of about 1.6. As the dotted lines in Fig.Eldemon- 
strate, the majority of the data sets follows the predicted 
scaling rather closely, whereas a the remaining ones show 
significantly slower decay. 

We have tried to analyze this finding further by looking 
at the distribution of squared displacements exhibited by 
all the particles, PMSB,t,{Sr'^), and its correlation with 
particle size. Here, i, is a fixed time, and the distribution 
is defined such that J PMSD,t^iSr'^)dr'^ = (5r^(i*). We 
have fixed t, such that (5r^ (i* ) = 1.25cP. The distribution 
^MSD develops a non-Gaussian peak centered around its 
average value, whose width increases upon increasing the 
packing fraction. In some cases, we did observe a two- 
peaked distribution at iy9 = 0.59, signifying that a certain 
amount of particles is displaced significantly less than 
average, i.e. that populations of 'fast' and 'slow' particles 
develop. This might be connected to the 'rare events' 
mentioned above, but we point out that this finding is 
unstable against improving the average over an increased 
number of simulation runs. 

From the a-scaling plot. Fig. 01 we infer the regime 
of a-relaxation dynamics. Note that for ip = 0.58, devi- 
ations from the a-master curve due to /3 relaxation are 
seen almost up to t « 100. Those will be analyzed later. 
Also shown in Fig. ^ are the MCT a-master functions. 
If evaluated at the same q as the simulation data, the 
description of the long-time dynamics is unsatisfactory, 
because the calculated stretching of the relaxation is too 
small. If we account for this error by shifting ^mct used 
in the calculations to higher values, we get good agree- 
ment, cf. the solid lines in Fig.0] Note that the deviations 
from the a-master curve set in at a time later than that 
where the ND and BD simulation results begin to over- 
lap: e.g., the q = 7.8 curve follows the a-master curve 
only for times t > 100, as can be inferred from Fig. 21 
Still, the BD and ND curves for that state collapse within 
our error bars already for t <; 20, cf. Fig. 13 This under- 
lines that the regime of structural relaxation identified in 
Fig. 121 is larger than that of the a-decay regime observed 
in Fig. 01 i.e. that both simulation data sets show some 
extent of the MCT /^-relaxation window. 

The values of ^mct used in the fits of Fig. 01 are 
gj^CT = 5.0 (9.13, 10.3, 15.13, 18.3, 20.87) for q = 4.0 
(7.8, 9.0, 13.8, 17.0, 19.8). These fit values are sug- 
gested by the analysis of the full curves pursued below, cf. 
Sec. IIVI Comparing the fitted wave-vector values ^mct 
to those of the simulations, deviations in q are in the 
range 10% to 17%, except for q — 4.0, where a 25% de- 
viation is needed to describe the a-master function. The 
way we have adjusted ^mct ensures that the stretching 
of the correlators is described correctly. In contrast, a 
fit of the plateau values with the a-master functions is 
difficult, since the latter are still not clearly visible in the 
simulation data even at (p = 0.59. This will become more 
apparent in Sec. IIVI 

In all cases, the fitted wave-vector values are larger 




FIG. 6: Example for Kohlrausch fits to the simulation data 
at (p — 0.58 (symbols; Brownian dynamics, dots: Newtonian 
dynamics), qd = 4.0, 7.8, 9.0, 13.8, 17.0, and 19.8. The fit 
range was t > 55. 



than the actual values, gMCT > q- Since the f^''^{q) 
giving the plateau values decrease monotonically from 
unity at g = to zero at q -^ oo, this fit result is 
equivalent to stating that the MCT-calculated plateau 
values appear too high. Furthermore, the half-width 
of the /*'^(g) distribution is an estimate for the in- 
verse localization length of a tagged particle. Hence the 
fit suggests that MCT underestimates the localization 
length of a tagged particle in the system slightly. There 
are two obvious reasons for such a mismatch in length 
scales: first, the softness of the particles in the simula- 
tion might, especially at high densities, give rise to some 
amount of particle overlap not possible in the HSS, ren- 
dering the effective localization of the particles slightly 
larger. According to Heyes [23, the soft-sphere system 
used in our simulations can be well described within the 
hard-sphere limit and an effective hard-sphere diameter 
4ff = /o°°(l - eM-PVc{r)])dr « d(l + 7e/36) w 1.016 
(where 7e ~ 0.577 is Euler's constant), which differs from 
d = 1 by less than 2%. But one has to keep in mind that 
the convergence of increasingly steep soft-sphere poten- 
tials towards the hard-sphere limit can be quite slow for 
the transport properties of the system 30). Second, a dif- 
ference in packing fractions between the simulation and 
the MCT calculation might become important in this re- 
spect. This arises, because within the PY approxima- 
tion for the DCF, the MCT master curve is evaluated at 
the corresponding value for the critical packing fraction, 
93^ « 0.516 < (^si^. As was pointed out in connection 
with Fig. ^ such a mismatch in ip will affect the aver- 
age particle distances, and thus an overall length scale. 
But since <Pmct ^ v'^ '-'^^ would expect this to lead to 
an overestimation of the critical localization length, con- 
trary to what we observe. 

Traditionally, stretched-exponential (Kohlrausch) 
laws. 



r(g,i)«^(9)exph(i/r(g))^('')] 



(9) 



are known to give a good empirical description of the a 
relaxation. Here, A{q) is an amplitude factor, r(q) the 





FIG. 7: Critical nonergodicity parameter f"''^{q) calculated 
from MCT for the one-component hard-sphere system with 
PY approximation (solid line), and for a five-component poly- 
disperse system (dashed line). Symbols are the amplitudes 
A{q) from Kohlrausch fits to the data, where error bars es- 
timate deviations depending on ip and BD/ND. The dash- 
dotted line indicates /"''^{q), but transformed with the wave- 
number shift applied in the discussion of the dynamical data 
(see text for details) . 



Kohlrausch time-scale of the a relaxation, and (3{q) < 1 
is called the stretching exponent. These parameters in 
general depend on the observable under study, and in 
particular on the wave vector q. Fig. |B1 demonstrates 
that the Kohlrausch laws can also be used to fit the a- 
relaxation part of our simulation data. The figure shows 
as an example the state ip — 0.58 for various wave vec- 
tors. One problem of the stretched-exponential analysis 
of the data is that the three parameters of the fit have 
a systematic dependence on the fit range. In particu- 
lar, one has to restrict the fitting to such large t that 
only a relaxation is fitted. For the fits shown, this range 
was chosen to be i > 55. The data deviate from the 
fitted Kohlrausch functions significantly only at shorter 
times; but there is a trend that these deviations set in 
just about the boundary of the fit range. This still holds 
if the fit is restricted to larger t only, and judging from 
the fit quality for the remaining relaxation alone, one 
cannot determine the optimal choice of the fit range. It 
is thus particularly difficult to extract the regime of a 
relaxation from the Kohlrausch fits alone. On the other 
hand, from the MCT fits shown in Fig. Qlwe expect cor- 
rections due to /3 relaxation to set in at about t « 100. 
This in principle gives an indication of the maximum fit 
range to chose. Yet, Kohlrausch fits restricted to t > 100 
did lead to an unsatisfactory scatter in the fit parameters 
A(q) and P{q). Thus an analysis of the a relaxation us- 
ing Kohlrausch fits will erroneously include parts of the 
f3 relaxation. 

This trend can be also identified comparing the ob- 
tained Kohlrausch amplitudes A{q) with the plateau val- 
ues f^''^{q). This comparison is shown in Fig. where 
the MCT results for f^'^{q) are included. They have 
again been determined using the PY approximation to 
5(q), but in agreement with Fig. [3 the values for f^''^{q) 
determined from MCT calculations based on the simu- 



FIG. 8: Stretching exponents f3{q) from fits to the BD simu- 
lation data at ip = 0.58 using Kohlrausch laws, Eq. (|5J. Error 
bars indicate deviations estimated from fits to ND data and to 
ip = 0.57. The dashed horizontal line indicates the b value cal- 
culated from MCT using the Percus-Yevick approximation for 
S{q), b ~ 0.583, and the dash-dotted line is b as determined 
from MCT with simulation-data input for S{q), b « 0.521. 



lated S{q) are indistinguishable from the ones shown on 
the scale of Fig. [7| For Kohlrausch fits to the a-master 
function, and more generally to the a-relaxation regime 
of the correlators only, A{q) < f'^'^{q) should hold. Re- 
calling the wave-number adjustment used in Fig. ^ we 
should even have A(g) < /''''^(gMCT(g))- This latter curve 
is included in Fig. [3 as the dash-dotted line, where the 
mapping q i~-> qmct was extended from the set of q an- 
alyzed in this text to all q via a quadratic interpolation. 
The relation A{q) < /*^^(g) is clearly violated for the fits 
here, especially at high q. It shows that the distinction 
between the a and /3 regimes from the simulation data is 
difficult; increasingly so with increasing wave number. 

It is reassuring that the Kohlrausch fits to BD and ND 
data yield values quite close to each other, apart from an 
overall shift in T{q). This holds, as long as the fit ranges 
are chosen such as to fit approximately the same part of 
the relaxation. In Fig. El the ND curves have been added, 
again shifted by a scaling factor in t given in the inset of 
Fig.O Note that, while in the ND curves one can identify 
a plateau from the data better than from the BD ones, 
still the KWW fits have a trend to give too high values of 
A{q). Thus one has to be careful when extracting plateau 
values from the simulation data by such an analysis, even 
if the data seem to give a clear indication of the plateau. 

The stretching exponents /3(q) from the KWW fits are 
shown in Fig. |S1 Again, we have included error bars in- 
dicating the deviations arising from fits to different ip or 
to ND as opposed to BD simulation data. /3(q) increases 
with decreasing q, and this increase is compatible with 
l3{q) — > 1 for (7 ^- 0, as expected from theory ^^. Ac- 
cording to MCT, P{q) should approach the von Schwei- 
dler exponent 6 as q ^ oo :32i. The value of b is cal- 
culated from the MCT exponent parameter A, and for 
the HSS using the PY-DCF is 6 w 0.583, shown as a 
dashed fine in Fig. |S1 We observe that the fitted /3(g) 
fall below this value for large q, even if the fits are less 
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FIG. 9: T{q) from Kohlrausch fits, Eq. @, to the BD data 
at (p = 0.58 (diamonds) and a,t (p = 0.59 (squares; scaled by a 
factor of 0.021). Error bars estimate the uncertainty from the 
fits, see text for details. The dotted line shows a 1/q^ law. 



reliable there, due to the low amplitudes A{q) at high 
q. To estimate the error of the theory prediction for b, 
we have also calculated this exponent from MCT using 
the simulated data as input for S{q). According to the 
values of A reported in connection with Fig. |21 we get 
b « 0.56 ± 0.04. The lower bound for b thus obtained is 
indicated in Fig. as the dash-dotted line. Taking into 
account this uncertainty, the behavior of the fitted f3{q) 
agrees well with what is expected from theory. In the 
further discussion, we will fix A to its value derived from 
the PY approximation, A — 0.735. Since the shape of the 
correlation functions in the /3 regime is in the asymptotic 
limit fixed by A, some deviations in the fits described 
below are to be expected in this time window. 



C. Analysis of Q-relaxation times 

The q-dependence of the a-relaxation times at fixed 
if can best be analyzed from the T(g) extracted from 
Kohlrausch fits. In Fig. |51 we report values for T{q) for 
such fits to the BD data at (p = 0.58 as the diamond sym- 
bols. If one instead fits the ND data, or data at ip = 0.585 
or 0.57, the q-dependence is the same up to a prefactor 
and up to small deviations. These deviations are indi- 
cated by the size of the error bars in Fig. |2| The data 
closely follow a l/q"^ dependence for small q, indicated by 
the dotted line. This is in agreement with earlier MCT 
predictions for the hard-sphere system [31|. For g — > cxd, 
one expects from MCT T(g) ~ q^^^''. But since b is close 
to 1/2, we cannot distinguish this behavior reliably from 
a 1/q^ law due to the noise of the data at large q. 

Fits to the BD data at ip — 0.59 reveal significant 
deviations from the behavior at (p < 0.585 at small q. 
This can be deduced from the square symbols in Fig. |51 
They have been scaled by a constant factor in order to 
match the value obtained from the (p = 0.58 fit at g = 7.8, 
since there the MCT analysis works best, as will be shown 
below. At smaller q, the increase of T{q) with decreasing 
q is suppressed for the ip = 0.59 data in comparison to 
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FIG. 10: Plots of TqT ''' for wave vectors qd = 4.0 (diamonds), 
7.8 (squares), 9.0 (up triangles), and 13.8 (down triangles); 
using 7 — 2.46. (a) Results for the simulation data using 
Brownian dynamics (open symbols) and Newtonian dynamics 
(filled symbols). The values of t for the latter have been 
multiplied by 4.5 for this plot. Dashed lines are linear fits to 
the Brownian dynamics data in the range 0.54 < ip < 0.58. 



the variation in T{q) observed for smaller p>. 

For a discussion of the density dependence of the a- 
relaxation time, the T{q) from the stretched-exponential 
fits are less reliable, since the Kohlrausch fits suffer from 
larger uncertainties at lower p, where the a- and /3- 
regimes are even less well separated. But we can op- 
erationally define a time scale Tq for the decay of the 
correlation functions as the point where the correlators 
have decayed to 10% of their initial value, </>*(<?, Tq.) = 0.1. 
For small enough q where f^''^{q) is still much larger than 
0.1, Ta is a useful measure for the a-process time scale. 
In the asymptotic regime, where a scaling holds, it fol- 
lows the a-scaling time defined within MCT up to a con- 
stant. Thus MCT predicts a power-law divergence of Ta 
close to p'^ of the form \p — p'^\~^ for not too large g, 
g < 15, say. To test this prediction, we plot in Fig. [TUI 

— 1/'7 . 

the quantity Tq , which should yield a straight line 
crossing zero at p'^ . Since the region of validity of this 
asymptotic result is not known a priori, a determination 
of the correct value of the exponent 7 on the basis of 
such rectification plot suffers from large errors. There- 
fore, let us fix 7 w 2.46, the value calculated by MCT 
using the PY approximation. Due to the uncertainty in 
determining A mentioned above, slightly different values 
of 7 cannot be ruled out. One gets 7 w 2.66 as an upper 
bound when using the upper bound for A given above 
for the MCT result based on the simulated SiqV This 
value of 7 is also quite close to what one gets JSlJ using 
the Verlet-Weis-corrected PY structure factor j3J]. Fig- 
ure El shows rectification plots for both Brownian and 
Newtonian dynamics simulation data for 4.0 < q < 13.8. 
For the latter, the values of Tq, have been multiplied by 
4.5, consistent with the shift in the overall time scale, cf. 
inset of Fig. |3| Fits to straight lines give p'^ values that 
are consistent with each other for q > 7.8 and both mi- 
croscopic dynamics, if one restricts the fit range to high 
enough p and omits the highest densities, where alpha 
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scaling breaks down, 0.54 < (p < 0.58 in our case. From 
this, one gets ip" « 0.594 ± 0.001. The data for q = 4.0 
give a somewhat higher value, ip'^ w 0.598 ± 0.001, again 
the same for Newtonian and Brownian dynamics. This 
difference is outside the error bars of the analysis and 
thus not in accord with MCT. Since the discrepancy is 
the same for both types of short-time dynamics, we con- 
clude that indeed the structural relaxation deviates from 
the MCT prediction systematically at small q. 

The range of distances e to the critical point, in which 
we can fit the time scales consistently with a power law, 
is roughly \ip — <p'^\lpf < 0.07. This agrees with what is 
expected from a discussion of the asymptotic MCT re- 
sults for the hard-sphere system 01 : For the time-scales 
extracted from the numerical MCT results, we have to 
restrict the linear fit to (/Jmct > 0.48, where the critical 
point is <(5^Qrp w 0.516. Below (pucT — 0.48, one finds 
deviations from the straight lines in the rectification plot; 
typically the results fall below the asymptotic straight 
line in such a plot. If one tries to fit a larger range in 
'/'MCT J the thus estimated p>'^ will be higher than the cor- 
rect one. For example, we get p'^ ~ 0.519 ± 0.0015 when 
fitting in the range (^mct > 0.4. It is reassuring that 
the deviations from the linear behavior seen in Fig. ^| 
for the simulation results occur in the same direction as 
found for the MCT results. 

At large q and at the highest packing fractions studied, 
the Tq from the simulations are systematically smaller 
than what is expected from the power-law extrapolation. 
This suggests that the local relaxation dynamics of the 
system very close to the transition would be faster than 
expected within the theory. However, the full theory 
analysis presented below suggests the opposite, indicat- 
ing that at these high q, the operational definition of r^ 
we have chosen for simplicity no longer works. 



D. /3-process analysis 

We now test some of the predictions MCT makes for 
the /9-relaxation regime, where the correlators remain 
close to their plateau values. The time window where 
the asymptotic solution holds, extends in t upon control 
parameters approaching the transition point, i.e. p ^ p'^ . 
The leading deviation from the plateau value is of order 
•\/|(t|, where a is called the distance parameter, and 



(7 = C • e, 



{p - p-)lp^ 



(10) 



in leading order is the linearized distance in control- 
parameter space. The leading-order asymptotic result 
is called factorization theorem, and it can be written for 
the tagged-particle correlator as 



<\>\q.t)=r'\q) + h\q)G{t) 



(11) 



Here, G(i) is a universal function depending only on the 
parameters A, a^ and a fixed 'microscopic' time scale to- 
The expansion is valid on a time scale ta — folo"!"^/'^"-* 




FIG. 11: (a) /3 analysis of the BD simulation data at </p = 
0.58: the curves marked by symbols show X(g, t) = (^^(g, t) — 
(t>%q,t'))/{(f,'{q,t') - 4>''{q,t")) with t' = 8.234, t" = 20.8075. 
Wave vectors are qd = 4.0 (diamonds), qd — 7.8 (squares), 
qd — 9.0 (up triangles), qd = 13.8 (down triangles), and qd = 
17.0 (circles). The dashed line is the equivalent of the MCT /3 
correlator, see text for details. The dash-dotted line indicates 
the plateau value estimated from the root of the (3 correlator, 
(b) Same for the ND simulation data, t' = 2.31 and t" = 
5.845. 



that diverges upon approaching the transition point. On 
this time scale, all wave-vector dependence is factorized 
off from the time dependence, and contained in the crit- 
ical amplitudes h''{q) and the plateau values /*'^((?). All 
parameters can be calculated within MCT, but as we 
have seen above, extracting them from the simulation 
data is not straightforward. Kob et al. |3J| have intro- 
duced a test of the factorization property that can be 
performed without fitting any of the g-dependent quan- 
tities to the data: they considered 



X{q,t) 



(t>'{q,t)~ct)%q,t') 
4>'{q,t')-<l>^{q,t") 



(12) 



If fixed times t' and t" are chosen inside the /3 regime, the 
factorization theorem gives X{q,t) = X{t) — xiG{t) — X2 
for t inside the (3 regime, with constants xi and X2 not 
depending on q. Therefore, if the factorization theorem 
holds, plotting X{q,t) for various q collapses the curves 
in the /3 window, without the need for fitting wave- vector 
dependent amplitudes and plateau values. 

We have performed this test for our simulation data 
for both BD and ND. Figure [TTk shows the results for 
the BD simulation at p — 0.58, with t' = 8.234 and 
t" — 20.8075. One observes that the data nicely collapse 
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for 5 < i < 40. Additionally, the figure shows the X{t) 
constructed from the MCT /3 correlator as a dashed line. 
Here, two constants xi and X2, and the time scale to- 
have been fitted. The value of A has been taken from the 
theory as explained above, A = 0.735. The same analysis 
is carried out for the ND data in Fig.lllb. Here, we have 



fixed t 



ND 



2.31 



t'/S.S and t'^j^ 



5.845 « i"/3.5. 



since a.t (p = 0.58 the shift in time scales between BD and 
ND is a factor of about 3.5, cf. inset of Fig.|3| Again the 
data collapse in an intermediate window 2 < i < 15. The 
upper end of this window is consistent with the one found 
for the BD analysis, i.e. 15 « 40/3.5. The lower end 
of the window where f3 scaling holds for the ND data is 
higher than what would correspond to the BD P window. 
Thus preasymptotic corrections are stronger in the ND 
case. The fit using the MCT /? correlator is not as good 
as it is for the BD case. Since the distance to the critical 
point docs not change between BD and ND, we have used 
the io- determined from the above fit to the BD data also 
here. Furthermore, since we have chosen t^j-, and ij^j-, 
in accordance with the values of t' and t" for the BD 
analysis, the constants xi and X2 should be the same; 
this is roughly fulfilled by our fit. 

The fits to both data sets have been performed such as 
to obey the 'ordering rule' for the corrections to (3 scaling 
13 : a curve that falls below another one for times smaller 
than the /3 window, will also do so for time larger than the 
f3 regime, since the corrections both at small and at large 
times are determined by the same q-dependent correction 
amplitudes. Thus the ordering of wave vectors on both 
sides of the scaling regime is preserved. This prediction 
of MCT is fulfilled by the BD simulation data, as can be 
seen in Fig. Illb . For the ND data, we cannot fulfill this 
ordering at both short and long times with reasonable 
fit parameters. At shorter times, one finds that e.g. the 
q — 7.8 and q — 9.0 curves rise above the /3-correlator 
curve, in violation of the ordering rule. We thus conclude 
that at this point, microscopic deviations set in for the 
ND simulations that are stronger than the ones in the 
BD data. The point qo where the corrections to /3 scaling 
change sign can be inferred from Fig. II II to be go ~ 9/d. 
It is independent on the type of short-time dynamics, 
in excellent agreement with the predicted universality of 
structural relaxation, and in particular the correction-to- 
scaling amplitudes. The numerical value of qo also agrees 
well with that found in an analysis of the MCT results 
for the tagged-particle correlator in a hard-sphere system 
[23, where the change occurs at qo,MCT ~ 9.3/d. 

The P correlator for short times approaches the critical 
power law, G{t) ~ i^". Comparing this asymptote with 
the fitted f3 correlator in Fig. ^] one finds that the t^°- 
law already deviates from the f3 correlator at i « 1 for 
the BD data, and at i « 0.3 for the ND data. Thus the 
critical decay cannot be identified from the simulation 
data. This is typical for most experimental data JJ . One 
thus has to be careful when extracting the exponent a 
from an analysis of the /3 relaxation. 

Let us from now on restrict the discussion to the BD 




FIG. 12: MCT (solid lines) and simulation results (symbols) 
for (j}"{q, t). For the simulation data, packing fractions are tp = 
0.50, 0.53, 0.55, 0.57, 0.58, 0.585, and 0.59 (from left to right), 
and qd = 7.8. For the MCT results, packing fractions have 
been adjusted to (^mct = 0.445, 0.47, 0.484, 0.499, 0.505, 
0.508, and 0.5135, and the wave number has been adjusted to 
QMCTd = 9.13; see text for details. 
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FIG. 13: MCT and simulation data for (j}^{q,t) with symbols 
and packing fractions as in Fig. 1121 but for qd — 9.0, adjusted 
in the MCT calculations to qmcTd = 10.3. 



data set. For the ND data, deviations from the MCT 
predictions occur in the early part of the P regime, and 
thus the theory can explain a larger part of the BD curves 
than it can do for the ND ones. For the analysis of the 
long-time universality of the dynamics outlined above, 
we conclude that these deviations are not a feature of 
the glassy dynamics. It is known that MCT treats the 
short-time dynamics insufficiently JT"! , and that Brownian 
dynamics typically stays closer to the MCT scenario for 
a larger time window than the corresponding Newtonian 
dynamics p7ll2^. 



IV. FULL MCT ANALYSIS 

We now turn to a full data analysis, i.e. a comparison 
of the complete simulation data with the solutions of the 
full MCT equations for a hard-sphere system. 

In the MCT picture, the glassy dynamics of the hard- 
sphere system is mainly driven by density fluctuations 
over the length scale of the mean nearest-neighbour dis- 
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FIG. 14: MCT and simulation data for (j)''{q.,t) with sym- 
bols and packing fractions as in Fig. 1121 but for qd — 13.8 

(ijMCTrf = 15.1). 
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FIG. 15: MCT and simulation data for cf)''{q,t) with sym- 
bols and packing fractions as in Fig. 1121 but for qd — 17.0 
(gMCxrf = 18.3). 



tance, i.e. with wave numbers close to that of the first 
sharp diffraction peak in S{q). We therefore begin the 
analysis by focusing on the data for q = 7.8. The re- 
sults of our fuU-MCT fits to the BD simulation data are 
shown in Fig. ^1 To achieve this and the following fits, 
we have adjusted (pucT for each curve and allowed the 
wave number qmct to vary slightly with respect to the 
correct value q. No other parameters have been adjusted. 
As noted above, the g-shift to some extent accounts for 
a mismatch in length scales between the simulation and 
the theory predictions. The adjustment of (pucT on the 
other hand accounts for the known error in (p'^. We will 
discuss the relation of the fitted ipuCT to the correct 
packing fraction p below. 

The fit shown in Fig. E| demonstrates that the the- 
ory can, with these modifications allowed, account for 
the dynamics of the hard-sphere system in a time win- 
dow of over four decades on a 10% error level. Only at 
larger times, t « 10^ in our units, i.e. at the highest pack- 
ing fraction studied, systematic deviations are observed. 
The simulation for this packing fraction shows slower dy- 
namics than expected from the theory. Also the shape of 
the final decay is different, as noted above in connection 
with a scaling. On the short-time side, the MCT de- 
scription works down to a time i « 1. At shorter times. 



it is still almost quantitative, but one observes a different 
curve shape. The simulation data appears more strongly 
damped than the MCT curves. This is to be expected, 
since neglecting the regular part of the memory kernel in 
Eq. (|3bl) will lead to such deviations at short times. We 
could have accounted for this partly by choosing a higher 
value of ly in Eq. Ij3b|l , but we have not done so since we 
are not concerned with the short-time dynamics in this 
work. 

Once the fit for q = 7.8 was completed to fix the 
empirical relation '^mct(</3), we have analyzed data for 
other wave numbers up to g = 30, hereby fixing the re- 
lation qMCTil)- Typical results for q < 17 are exhibited 
in Figs. El to ^1 Note that the only parameter that 
was adjusted for these fits is the wave number, quCT- 
This way. Figs. 1131 to 1151 demonstrate how MCT is able 
to reproduce the wave-vector dependent changes in the 
structural-relaxation window. At even higher q, it be- 
comes too difficult to judge the fit quality, since the f{q) 
are close to zero there. Connected with this is the grow- 
ing infiuence of the microscopic regime, t < 1, on the 
main part of the decay of the correlators with increasing 
q. For g = 17 and the highest packing fractions, already 
about 60% of the decay of 0^(q, t) from unity to zero are 
made up of such microscopic dynamics. Consequently, 
the errors made in its description are to be seen more 
exphcitly in Fig. [THl than in Fig.lT^ 

Apart from this, also the fits at g > 7.8 using the fuU- 
MCT results are quite satisfactory in the time window 
1 < t < 10^, save the highest simulated density, for which 
errors are largest and extend down to t w 10 for q = 13.8 
and 17, cf. Figs. 1141 andllSI One notices a trend that the 
ck-relaxation dynamics becomes slower in the simulation 
data than expected from the MCT fits, i.e. the local dy- 
namics is slower than one estimates in the theory. The 
finding can probably not fully be attributed to the incor- 
rect structure-factor input used, since a recent study of 
binary hard-sphere mixtures reported a similar discrep- 
ancy for the g-dependence of the a-relaxation times even 
when basing MCT on the 'correct' simulated S{q) as in- 
put ;8i] . The same trend is also apparent in the fuU-MCT 
analysis of a binary Lennard- Jones mixture |3 . 

Having established the overall quality of the MCT de- 
scription for the structural dynamics on length scales 
smaller and comparable to the typical particle-particle 
distance, let us now investigate the adjustment in ipMCT 
needed to achieve this level of agreement. A plot of i^mct 
vs. p is shown in Fig. Iltil fdiamond symbols). The figure 
reveals that the relation is close to linear, and thus the 
nontrivial variation of the relaxation curves close to the 
singularity p"^ has not been put in 'by hand' through the 
fitting parameter pmct- We can estimate the correct 
value of the glass-transition packing fraction by a lin- 
ear fit to the <y9MCT-versus-(^ curve. Using the calculated 
value V'mct ~ 0.516, we get p'^ « 0.594. This value is 
nicely consistent with the one obtained from the a-scale 
analysis of the data, cf. Fig. (TUI and also from an earlier 
analysis of the simulations ,3^] . Note that it differs from 
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FIG. 16: Plot of (fiMCT vs. ip for the fits shown m Figs. 1121 - 
1151 1171 and 1181 (diamonds). The dashed line is a linear fit, 
ipMCT ~ O.Slip + 0.037. The circles are for the independent 
fit of the MSD, Fig. ^| The dotted horizontal line indicates 
the calculated critical point, if'^ ~ 0.516. 
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FIG. 18: Comparison of the mean-squared displacements 
5r^{t) from simulation and MCT calculations; all fit param- 
eters have been taken from Fig. 1121 In addition, the MCT 
curves have been multiplied by 1.1 to account for an error in 
localization length; see text for details. 




FIG. 17: MCT and simulation data for (f)''{q,t) with symbols 
and packing fractions as in Fig. 1121 but for qd = 4.0 [qmcTd — 
5.0). 



the result obtained from MCT based on the simulated 
S{q), (f'^ « 0.585 by less than 2%. The slope of the linear 
fit in Fig. 1161 is not equal to unity, and its zero is shifted. 
If one considers the connection of the distance parame- 
ter a of MCT to the control-parameter distance e, this 
translates into an error of the leading-order constant of 
proportionality C in Eq. (|10|) . From Fig. ^| we conclude 
that the value calculated from MCT, Cmct ~ 1-54 |2]j is 
in error by about 20%, C « 1.2. 

For small q, the MCT description of the data shows 
larger quantitative errors, while it remains qualitatively 
correct. This is exhibited by the fits done for q — 4.0, 
Fig. El Again, only ^mct was allowed to vary, while 
the (pMCT have been determined from the q = 7.8 fit 
shown in Fig.[Tl While for this latter fit, the MCT fits 
reproduce the a-relaxation times rather well, this is not 
the case for the q — 4.0 fit at (/5 > 0.57. Instead, one 
observes a systematic trend for the simulation data to 
decay increasingly faster than the MCT curves with in- 
creasing if. In addition, the deviations observed in the 
/3-relaxation window, while still within a 10% level, are 
larger for q = 4.0 than they are for q > 7.8. Even more, 
it was necessary to allow for a 25% deviations between 



Qmct and q, whereas this deviation was less than 17% for 
all other fits. This last finding suggests that the /^'^(g)- 
versus-q curve calculated within MCT is too broad. It 
is common to express deviations of the (f>^''^{q,t)-veis\is-q 
curve at fixed t from a Gaussian at small q in terms of 
the non-Gaussian parameter, a2{t). These non-Gaussian 
corrections are reproduced qualitatively correct by MCT, 
but with an error in magnitude. One finds a2{t) < 0.3 
in the theory [23 , while for our simulation, a2 (t) reaches 
values up to 2.5 in both BD and ND, which is in agree- 
ment with similar simulation results for other systems 
[33 ■ But note that for times where 0''(g,t) is close to its 
plateau value /"'^(g), both the MCT and the simulation 
value of a2{t) are positive. Thus the underestimation of 
a2 in the theory would let the /''•°(g)-versus-g curve ap- 
pear too narrow, opposite to what is observed from our 
fits. We thus conclude that non-Gaussian corrections as 
expressed through 02 (i) and the quantitative error MCT 
makes in expressing them cannot be alluded to to ex- 
plain the deviations observed at g = 4.0. Let us point 
out that the deviations discussed above do not depend 
significantly on the fact that we have based the MCT 
calculation on the PY-DCF. Using the simulated struc- 
ture factor with MCT gives correlation functions (t)^{q,t) 
that behave qualitatively as the ones shown here. 

It is instructive to extend this analysis towards the 
mean-squared displacement data. Since the MSD is given 
through a memory kernel that basically is a g — > limit 
of the tagged-particle-correlator memory kernel, Eq. iQ, 
its analysis can be viewed as the g ^ extension of the 
above fitting procedure. 

We report the BD simulation data for the MSD, to- 
gether with the MCT curves according to the correlation 
functions shown above, in Fig. ^1 For the MSD, no 
wave number is to be adjusted, and in this sense, the 
MCT results of Fig. E| are not fitting results, but rather 
consequences of the analysis done for q — 7.8, shown in 
Fig- El We have, however, adjusted a global length scale 
in this plot, in order to better fit the localization plateau 
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FIG. 19: Comparison of the mean-squared displacements 
Sr^{t) from simulation and MCT calculations. In this plot, 
values <^MCT used for the MCT calculations have been ad- 
justed to fit the long-time diffusion regime of the simulated 
Sr^{t); the values are <^mct = 0.438, 0.46, 0.478, 0.494, 0.502, 
0.503, and 0.511 for <^ = 0.50, 0.53, 0.55, 0.57, 0.58, 0.585, 
and 0.59, respectively. 




FIG. 20: Self-diffusion coefficients D as a function of the 
packing fraction (p, as obtained from the Brownian dynam- 
ics simulation (diamonds, with connecting lines to guide the 
eye) and from MCT calculations (solid line) . The dashed line 
indicates the MCT asymptote, D oc {ip" — f)'^ ■ 



visible in the data. The MCT curves have been scaled 
up by a factor of 1.1, which accounts for a 5% underes- 
timation of the localization length by the theory. Note 
that at short times, i < 1, the description of the data 
using MCT is of similar quality as discussed above. In 
particular, the MSD plot reveals that the BD simulation 
still resembles a Newtonian short-time dynamics, though 
strongly damped: the MSD roughly follows a 5r^ ^ P 
law for 10^^ < t < 10^^, and not a Sr^ ^ t law as would 
be expected for short-time diffusion in a strictly Brown- 
ian system. Theory and simulation do not match at short 
times because of the scale factor applied. For long times, 
a qualitatively similar picture emerges as for q = 4.0, 
regarding the variation of the a-relaxation times with 
If, now showing through a displacement of the long-time 
diffusive straight lines in the 6r^{t) plot. As for q = 4.0, 
the quality of the MCT description of the /3-relaxation 
regime similarly is worse than for q > 7.8. But for the 
MSD, all deviations are larger than for q — 4.0, especially 
in the a-relaxation regime. 

Before investigating the a-relaxation regime in more 
detail, let us note that, all deviations taken aside, the 
shapes of the MSD curves as predicted by MCT are qual- 
itatively correct. To substantiate this statement, Fig. 1191 
shows an independent fit of the MSD using MCT. In- 
stead of fixing the (/SMCT-versus-tp relation from the data 
at q = 7.8, as was done above, in this case this relation 
was determined from the MSD alone. It is noteworthy 
that by correcting the error in the a-relaxation time scale 
observed before, also the description of the /3-relaxation 
window improves. In particular, we did not scale the 
MCT results as we have done in Fig. ^1 The values 
'PMGTif) used in Fig. El are reported in Fig. El as the 
circle symbols. They also lie on a straight line, which is 
shifted downward somewhat with respect to the original 
relation deduced from the above fits. In turn, an esti- 
mation of ip'^ from the mean-squared displacement, i.e. 
from the diffusivities, yields a value that is too high, viz.: 



(p'^ « 0.598 ± 0.003. This is a typical result also found 
in other simulations 8], but not in accord with MCT. 
But note that the estimations for ip"^ are quite close to 
each other, so the deviations can be regarded as indi- 
cations of error margins. As a side remark, let us note 
that the parameters from the independent fit to the MSD 
presented in Fig. ll9l could be used to improve the descrip- 
tion of the q = 4.0 case somewhat, but not completely. 
One concludes that the MCT description of the single- 
particle structural relaxation smoothly deteriorates for 
decreasing q. 

The deviations in the a-relaxation regime, i.e. the long- 
time diffusive regime, that arise for the MSD can be 
analyzed more clearly by looking at the long-time self- 
diffusion coefficient D{(p) itself. Here, D has been de- 
termined from the simulations by the Einstein relation, 
Sr^ ~ 6Dt for large t, at times t where 6r^ is of the order 
of 10 squared particle radii. The results for the BD sim- 
ulations are shown in Fig. 1201 as the diamond symbols. 
In contrast, the MCT calculation, plotted in the figure 
as a solid line, shifted according to the relation ipucTiv) 
used in the above discussion, systematically falls below 
these values. The relative error in D is less than 20% up 
to (fi Ri 0.55, increases to 80% at (p = 0.58, and reaches a 
factor of 4 at ip — 0.59. The two curves could be matched 
within the error bars by further shifting the MCT results 
along the ip axis by less than 1%, which is basically what 
has been done in Fig. El But let us stress that there is 
no theoretical justification for doing so. 

MCT predicts a power-law asymptote for D with the 
same exponent 7 that applies for the divergence of the 
a-relaxation times, D ex \(j\'^ ■ This asymptote is included 
in Fig. 1201 for the MCT results as a dashed line. It is also 
possible to fit the simulation data with such a power-law. 
We have restricted these fits to (p > 0.55 and have omit- 
ted the value at (/? = 0.585. If we fix the exponent to the 
theoretical value, 7 w 2.46, we get reasonable agreement 
with a fitted (p'^ « 0.599, in agreement with the above 
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FIG. 21: Product D -T^q) at wave vector qd = 7.8 for the BD 
simulation (crosses, connected by lines to guide the eye) , and 
for the MCT curves (solid line). The latter curve has been 
transformed along the horizontal axis according to Fig. Ilfcil 
The inset shows a magnification of the MCT result versus 

</5MCT. 




FIG. 22: Averaged static structure factor S{q) for a monodis- 
perse system (dashed lines) and the polydisperse system stud- 
ied in this work (solid lines) at the same state point. The 
main figure shows the simulation results at i/p = 0.54. The 
inset shows the Percus-Yevick S{q) at the respective MCT- 
critical packing fractions ^Smcti for both the monodisperse 
(dashed) and a five-component system (solid lines). 



observation. If we, on the other hand, also determine 
7 from the fit to D{(p), the result is cp'^ « 0.597 with 
7 « 2.24. It is a typical observation from simulation and 
experiment is that an independent determination of 7 
from the diffusion coefficient yields a different value than 
that obtained from the analysis of the density correla- 
tors I33I • From the comparison of the MCT results with 
the asymptotic prediction in Fig. 1201 it is, however, clear 
that the asymptotic law only holds for D < 10~^, i.e. for 
(f > 0.58 for our simulations. Thus a large part of the fit- 
ted simulation data is outside the asymptotic regime for 
this power law, and the fit yields an effective exponent 
rather than the true 7. 

The above results indicate that with (p increasing close 
to v?'^, a decoupling of the diffusion time scale, as seen 
from the mean-squared displacement, from the density- 
fluctuation-relaxation time scale, as seen in the density 
correlation functions, sets in. This can be illustrated 
without referring to any fits by plotting the product 
Dr^q) of the diffusion coefficient D and the a-relaxation 
time scale |33j|. For the latter, let us choose the value 
obtained for q = 7.8, as a typical representative of the 
local-order length scale. Figure ETl shows as symbols the 
results from the BD simulation. One notices an increase 
in Dt by a factor of 2 to 3 within the density range cov- 
ered by the simulations. We have also checked that this 
holds similarly for q = 4.0 and q = 9.0. The correspond- 
ing MCT result is shown as the solid line in Fig. [21 which 
is magnified in the inset of the figure. Here, the product 
Dt also increases with increasing ip close to t^'^, but only 
on the order of 10%. One thus concludes that there is a 
rather large quantitative error in this quantity, although 
not necessarily a qualitative one. MCT predicts that Dt 
approaches a finite value as (^ — > i/j"^. As to whether Dt 
diverges or stays finite at Lp'^ in the simulation, our data 
remain inconclusive. Note that the values for the highest 
two packing fractions simulated are relatively uncertain, 
as the scatter in Dt indicates. 



V. POLYDISPERSITY EFFECTS 

Up to now, we have neglected the fact that the sim- 
ulated system is not strictly a single-component system. 
Instead, some size polydispersity was needed in order to 
avoid crystallization. In this section, we give a brief ac- 
count of how much we expect this small polydispersity 
to affect the results discussed above. 

We first examine the infiuence of polydispersity on the 
equilibrium fluid structure. To this end, we have sim- 
ulated a monodisperse system of the same soft spheres 
as used in the polydisperse simulations. We found such 
simulations possible for packing fractions up to (^ sa 0.54, 
above which crystallization as monitored by Qe sets in 
rather quickly. The static structure factor S{q) at this 
state point is compared to the one from the polydisperse 
system at the same density in Fig. |221 As expected, the 
polydisperse system exhibits less pronounced ordering, 
visible in reduced oscillation amplitudes in S{q). The 
effect is well explained by the PY approximation, as 
the inset of Fig. |221 demonstrates. There, the (total- 
density) structure factor for the one-component hard- 
sphere system is compared with that obtained from the 
five-component mixture introducted in Sec. IIIBI One 
might expect the visible differences in the monodisperse 
and the polydisperse S{q), however small, to affect the 
MCT results for tp'^. This would be true if both systems 
were treated as one-component systems. But it is not 
necessarily true for a full calculation of multi-component 
MCT, using the full matrix of partial structure factors 
instead of only the averaged S{q), as we will discuss now. 

Let us compare the results obtained for /"•'^{q) for the 
one-component system with those of the five-component 
system mentioned above. This comparison is included in 
Fig. 13 where the dashed line shows the averaged fp^{q) 
according to Eq. ||HJ). The thus obtained curve is slightly 
narrower than the solid curve, representing the result of 
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the one-component calculation. Accordingly, the average 
localization length increases slightly, by about 4%. From 
the above discussion we conclude that this is a change 
in the right direction, but not enough to account for the 
wave-number shift we had to apply to describe the den- 
sity correlation functions with the one-component sys- 
tem. Comparing with the Kohlrausch amplitudes A{q) 
also shown in Fig. {7\ we note that the intrinsic error in 
determining the plateau values from the data is larger 
than the differences between the two MCT curves. 

The values of ip'^ obtained from the MCT calculations 
with one, three, and five components show only minor dif- 
ferences. While the one-component result is ip'^ « 0.5159, 
we get <Pg^3 « 0.5153, and (fi^^^ ~ 0.5154. Similarly, the 
exponent parameter only changes slightly between these 
systems: from A « 0.735 in the one-component system to 
As=5 ~ 0.737 for the five-component case. These changes 
in if'' and A are significantly smaller than the uncertainty 
in these quantities coming from the approximation used 
for the static structure factor. Note that the value of (/3^ 
decreases slightly in the multi-component systems men- 
tioned. This is consistent with recent MCT predictions 
for a two-component system |2(tJ . where it was found 
that for size ratios dsmaii/'^iargc ^ 0.8 the critical point 
If'' slightly decreases compared with the one-component 
system. Only when the size ratio became more extreme, 
c^smaii /enlarge < 0.6, Say, did the MCT calculations show 
a notable effect on ip'^. In this latter case, the values 
of if'' were found to be larger in the mixture than in 
the one-component system. This increase is commonly 
expected for polydisperse systems. But from our cal- 
culations we conclude that such polydispersity-induced 
fluidization does not occur for the present polydisperse 
size distribution, which in particular lacks any large- or 
small-radius 'tails'. 

A further comparison to the predictions of the multi- 
component MCT can be made by binning the particles 
of the simulation according to their size into a differ- 
ent number of bins. Let us demonstrate this for the 
case of three bins, a — small (radii 0.9 < dsmaii < 
0.96667), medium (0.96667 < dmcdium < 1.03333), and 
large (1.03333 < diargc < 1-1). The thus obtained three 
tagged-particle correlation functions (p^ {q, t) can be com- 
pared to the three tagged-particle correlation functions 
amenable to the MCT calculation in the three-component 
system. Figure |23K shows as symbols the results from a 
three-bin analysis of the BD simulation data aX ip — 0.58 
and q = 7.8. One notices that the largest particles 
show the slowest decay, while the smallest particles de- 
cay fastest, as is intuitively expected. In the /3-relaxation 
window, one finds an ordering of the plateau values from 
small to large particles, where the smallest particles show 
the smallest plateau. Again this follows the expectation 
that the particles are localized more tightly the larger 
they are. These qualitative features are in agreement 
with the MCT results for the three-component mixture, 
as can be inferred from the symbols Fig. 123b . The order- 
ing of the plateau values is indicated by the horizontal 




FIG. 23: (a) Correlation functions (j>a{<l,t) for the BD sim- 
ulation, binned into three particle sizes, a — small (dia- 
monds), medium (plus symbols), and large (circles); see text 
for details. The data refer to packing fraction ip = 0.58 and 
qd — 7.8. The averaged quantity (/>pd('7, t) analyzed in the dis- 
cussion in detail is plotted as a solid line but coincides with 
the a — medium curve on the scale of the figure. The solid 
and dashed lines decaying at shorter times are the results for 
</5 = 0.54 from the simulation of the polydisperse and the 
monodisperse system, respectively, (b) As in (a), but results 
from the MCT calculations for a three-component mixture at 
packing fraction ip = 0.505 and qd = 7.8. Again, the aver- 
aged <^pd(9,i) is included as a solid line and is hidden by the 
a = medium curve. The solid line decaying at shorter times 
is the averaged result from the three-component mixture at 
ip = 0.4. For comparison, one-component results ai p — 0.42 
and p = 0.507 are included as dashed lines, the latter being 
obscured by the polydisperse-averaged p — 0.505 result. 



solid lines which represent the results for fa''^{q). Note 
that the differences in the relaxation curves for the three 
components are more pronounced than in the binned 
analysis of the polydisperse simulation, which might be 
related to the fact that the particle size distribution in 
the simulation is continuous. In the MCT calculation, the 
alpha-relaxation time of the large particles, measured by 
0[* j,((7, r) = 0.1, is at the wave number chosen slower 
by a factor of about 1.48 compared to that of the small 
particles. This change is slightly higher than in the sim- 
ulation, where the same trend applies with a factor of 
about 1.25. 

The unbinned correlation function from the ip — 0.58 
simulation, averaged over all particles, is included for 
in Fig. 123b . but on the scale of the plot, it coincides 
with the correlation function for the a — medium bin. 
Thus, in a sense, polydispersity effects 'average out' in 
this quantity. A similar effect applies for the MCT re- 
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suits, but here, the small difference in the critical packing 
fraction induced by the polydispersity leads to a shift 
in the relaxation time scales close to the glass transi- 
tion. Still, the one-component correlator calculated at a 
slightly higher packing fraction, </? = 0.507, agrees with 
the a — medium correlator from the three-component 
mixture at 1^9 = 0.505 on the scale of the figure. This 
agreement is not too surprising, because the MCT pa- 
rameters quantifying the slow relaxation features apart 
from if'^ do not change significantly between the monodis- 
perse and the three-component system, as was mentioned 
above. At lower ip, however, small differences become 
more apparent. This can be seen in the three-component 
ip = OA correlator shown in Fig. 123b as a solid line. It 
compares well with the result from a monodisperse cal- 
culation (the dashed line in Fig. I23b'l. but at (p = 0.42; 
and one notices somewhat different curve shapes. Again, 
these polydispersity effects are even smaller in the simu- 
lations. The solid line in Fig. 123b shows the simulation 
result for the polydisperse system at ip — 0.54, which is 
compared to the result from the monodisperse simulation 
at the same density, shown as a dashed line. Here, the 
agreement between the two systems is even better; and 
note that we did not have to adjust the packing fractions 
in this comparison. 

Thus it appears that this way of representing the poly- 
disperse system as a three-component mixture leads to a 
systematic overestimation of polydispersity effects. For 
the binary mixtures studied in Ref. iq, it was found 
that MCT even underestimates the size of the observed 
mixing effects. If this applies also to our case, the 
over-estimation of polydispersity effects by the three- 
component approach would be even stronger. 



VI. CONCLUSION 

We have performed Newtonian (ND) and strongly 
damped Newtonian dynamics (BD) simulations of a poly- 
disperse quasi-hard-sphere system close to the glass tran- 
sition. The wave-vector dependent tagged-particle cor- 
relation functions and the mean-squared displacement 
curves have been analyzed using the mode-coupling the- 
ory of ideal glass transitions (MCT), in order to provide a 
stringent test of the complete theory for a reference case. 

To ensure that the simulation data show all the quali- 
tative features predicted by MCT close to the glass tran- 
sition, we have analyzed both the ND and BD data in 
terms of a- and /3-scaling, cf. Figs. 0] and ^] This al- 
lows us to identify the time domain, where one can ex- 
pect MCT to give a quantitative description of the data, 
t > 10 in our units. In particular, both ND and BD agree 
at long times up to a trivial time scale. This universality 
of the structural relaxation is predicted by the theory, 
and fulfilled in great detail by our simulation data. In 
particular, the /3-scaling parameters and those qualita- 
tive features of the correction-to-scaling amplitudes we 
could test, are independent on the type of short-time dy- 



namics. Similarly, an analysis of the a relaxation with 
stretched-exponential fits demonstrates that the wave- 
number dependent shape of this relaxation is in agree- 
ment with what one expects from MCT. Other param- 
eters, as for example the critical-decay power law pre- 
dicted as an asymptotic MCT feature, cannot be ex- 
tracted from the simulation data. The analysis reveals 
that the highest density studied in our work, ip = 0.59, 
shows systematic deviations from the MCT predictions, 
and can thus not be explained by the theory. 

The main purpose of this paper is to compare the sim- 
ulation data to the full solution of the MCT equations. 
Leaving aside the small difference between the simulated 
polydisperse soft spheres and a true hard sphere system, 
this comparison is, in principle, free from any parame- 
ters. Since MCT is an approximate theory, one expects, 
however, deviations that can be accounted for by allow- 
ing some of the parameters to vary. 

We are able to achieve very good agreement between 
theory and simulations if we allow for a smooth mapping 
of packing fractions, ip 1— > ipmCT, and a similar mapping 
of wave numbers, q i~-> (/mct- The reasons for the needed 
adjustments are well understood. First, the critical point 
for the (ideal) glass transition calculated within is too 
low. This is compensated by mapping (p. The mapping 
turns out to be almost linear, and hence inessential in 
order to understand the slow relaxation close to the glass 
transition as a function of the distance to this transition. 
Second, we observe a small mismatch in length scales be- 
tween the simulation and the MCT results. This can be 
accounted for by mapping q. The difference in length 
scales is typically of the order of 15%, and only about 
5% for the localization length estimated from the mean- 
squared displacement. It is possible that these discrep- 
ancies are to some extent due to the slight softness of the 
simulated system. 

Given these parameter mappings, MCT is able to de- 
scribe the BD-simulation data over most of the density 
range studied quantitatively on a 15% level, as demon- 
strated in Figs. E] to ^] This extends down even to 
relatively short times, i « 1, and over a large range 
of length scales, from the nearest-neighbour distance 
down to qd « 20. At larger length scales (smaller 
wave numbers), stronger deviations set in, which are 
most pronounced in the long-time diffusive regime of the 
mean-squared displacement, but also in the /3-relaxation 
regime, cf. Figs. 1171 andll8l 

One has to keep in mind that the kind of comparison we 
have performed here is influenced by three conceptually 
distinct error sources: (i) deviations due to the compar- 
ison of a slightly polydisperse system in the simulations 
to a strictly monodisperse one in the theory; (ii) devi- 
ations due to incorrect structure-factor input to MCT; 
and (iii) deviations inherent to the MCT approximation. 
In order to shed more light on the quality of the MCT 
approximation itself, we have tried to disentangle these 
three error sources. 

The influence of polydispersity in the studied system is 
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negligible, as we have pointed out in Sec. by compar- 
ing to a three-component and a five-component system 
mimicking the polydispersity distribution imposed in the 
simulations. 

On the other hand, the second error source, due to ap- 
proximations made in describing the static equilibrium 
structure, has to be considered carefully. We have cho- 
sen to base most of our discussion on MCT results cal- 
culated from the Percus-Yevick structure factor for the 
hard-sphere system, because this is closest to a first- 
principles calculation. However, if one bases MCT on 
the structure factor obtained from the simulation, one 
can improve the description of the data. Most promi- 
nently, this infiuences the prediction of the critical point, 
which shifts from (^mct ~ 0.516 to (/^mct ~ 0.585, and 
thus surprisingly close to the experimentally determined 
value. At the same time, many of the predictions based 
on the PY structure factor remain quantitatively true, 
such as the shape of the /'^(g)-versus-g distribution, or 
the asymptotic shape of the correlation functions. This 
finding to some extent justifies our approach of adjusting 
the packing fraction ip in the PY-based MCT calcula- 
tions. In principle, a further error source connected with 
the static-structure input is the factorization of three- 
point static correlations in the MCT vertices, Eqs. (|3d|) 
and H5c|l . But we expect this purely technical approxima- 
tion to have small influence for our system, as simulation 
studies of a binary Lennard- Jones mixture |5j suggest for 
systems dominated by hard-core repulsion. 

The remaining discrepancies between the simulation 
results and the MCT predictions for the hard-sphere sys- 
tem are likely to be the ones giving information about 
the quality of the MCT approximation itself. These are: 

(i) The wave-vector variation of relaxation times. This 
is less pronounced in MCT than it is in the simulations. 
For large wave numbers, the BD simulation shows slower 
relaxation than expected from the theory, while at small 
q, the relaxation is faster than predicted by MCT. This 
indicates an error of the theory in capturing the length- 
scale dependence of the dynamics. The error at small q 
might be more severe, and is most dramatic when one 
considers the mean-squared displacement, i.e. the diffu- 
sion coefficient. The theory predicts that all structural 
relaxation time scales are coupled close to the glass tran- 
sition. This implies that the product of the diffusion 
coefficient and a typical intermediate-length-scale relax- 
ation time, D ■ r, should approach a constant when (p 
approaches ip'^. For finite pf — p, MCT predicts a smooth 
variation that is in qualitative agreement with the sim- 
ulation results, cf. Fig. [2] But the magnitude of this 



variation is underestimated by a factor of 2.5. This can 
be viewed as a quantitative error that has, however, a 
large impact on the description of D or the relaxation 
times at small q. We could not test whether the simula- 
tion behaves qualitatively different to the MCT predic- 
tion as (/9 — > (ys'^ due to obvious constraints. In general, 
our results show that the improper treatment of time- 
scale decoupling within MCT is not peculiar to the diffu- 
sion coefficient itself, but rather sets in smoothly at small 
q in the (j)''{q,t). 

(ii) At short times the description of the relaxation 
curves with MCT is insufficient. This is known since 
long, but it remains an important question at how large 
times deviations can still be seen. In our comparison of 
strongly damped Newtonian dynamics, is appears that 
the tagged-particle correlators can be fitted quite well 
even down to i ~ 1, i.e. almost 'microscopic' time scales. 
But a comparison with undamped Newtonian dynam- 
ics simulations reveals that there, short-time corrections 
can occur even for relatively large times, up to t w 10 
in our case. This can provide an explanation for recent 
observations stating that for a binary mixture obeying 
strongly damped colloidal dynamics, the MCT descrip- 
tion extended quantitatively down to surprisingly short 
times 9] , whereas a similar comparison of Newtonian MD 
data was satisfactory only in the a-relaxation regime |g] ■ 

(iii) At the highest packing fraction analyzed in the 
present work, more dramatic deviations between the sim- 
ulation results and MCT occur. They are most easily 
observed as a violation of a scaling, and a different scale 
behavior of the corresponding relaxation time. Our simu- 
lations hint towards possibly rare events that induce this 
behavior. But we have not been able to establish this 
within reasonable statistics. 
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